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Abstract 

Two-phase  flow  and  transport  of  reactants  and  products  in  the  air  cathode  of  proton  exchange  membrane  (PEM)  fuel  cells  is  studied 
analytically  and  numerically.  Single-  and  two-phase  regimes  of  water  distribution  and  transport  are  classified  by  a  threshold  current  density 
corresponding  to  first  appearance  of  liquid  water  at  the  membrane/cathode  interface.  When  the  cell  operates  above  the  threshold  current 
density,  liquid  water  appears  and  a  two-phase  zone  forms  within  the  porous  cathode.  A  two-phase,  multicomponent  mixture  model  in 
conjunction  with  a  finite-volume-based  computational  fluid  dynamics  (CFD)  technique  is  applied  to  simulate  the  cathode  operation  in  this 
regime.  The  model  is  able  to  handle  the  situation  where  a  single-phase  region  co-exists  with  a  two-phase  zone  in  the  air  cathode.  For  the 
first  time,  the  polarization  curve  as  well  as  water  and  oxygen  concentration  distributions  encompassing  both  single-  and  two-phase  regimes 
of  the  air  cathode  are  presented.  Capillary  action  is  found  to  be  the  dominant  mechanism  for  water  transport  inside  the  two-phase  zone  of 
the  hydrophilic  structure.  The  liquid  water  saturation  within  the  cathode  is  predicted  to  reach  6.3%  at  1.4Acm~2  for  dry  inlet  air. 
©  2001  Elsevier  Science  B.V.  All  rights  reserved. 
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1.  Introduction 

Water  management  is  critical  for  achieving  high  perfor¬ 
mance  of  proton  exchange  membrane  (PEM)  fuel  cells.  As  a 
vessel  of  protons,  the  polymer  electrolyte  membrane 
requires  sufficient  water  to  exhibit  a  high  ionic  conductivity. 
During  fuel  cell  operation,  water  molecules  migrate  through 
the  membrane  under  electro-osmotic  drag,  fluid  convection, 
and  molecular  diffusion,  making  it  difficult  to  retain  a  high 
water  content  within  the  membrane.  Generally,  humidifica¬ 
tion  is  applied  to  the  inlets  of  the  anode  and/or  cathode  in 
order  to  supply  water  to  the  membrane  region.  On  the  other 
hand,  water  is  generated  at  the  cathode/membrane  interface 
due  to  the  electrochemical  reaction  of  H+/02-  If  the  water 
generated  is  not  removed  from  the  cathode  at  a  sufficient 
rate,  cathode  flooding  may  result  and  the  oxygen  gas  trans¬ 
port  is  hindered.  Thus,  a  relatively  dry  air  at  the  cathode  inlet 
is  sometimes  helpful  to  remove  excessive  water,  in  parti¬ 
cular,  for  the  air  cathode  of  direct  methanol  PEM  fuel  cells. 
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Modeling  water  transport  and  distribution  in  PEM  fuel 
cells  has  been  attempted  in  recent  years.  Bernardi  [1] 
analyzed  the  effects  of  different  operating  conditions  on 
the  water  balance  in  fuel  cells.  Wang  and  Savinell  [2] 
included  the  electrode  structure  effects  on  water  transport 
in  the  catalyst  layer.  With  a  pseudo  two-dimensional  mem¬ 
brane-electrode  model,  Fuller  and  Newman  [3]  simulated  a 
PEM  fuel  cell  operation  under  the  condition  of  moist  gas  in 
electrodes.  They  obtained  the  water  distribution  in  the 
membrane  and  the  net  water  flux  across  it.  Nguyen  and 
White  [4]  compared  three  humidification  designs  using  a 
combined  heat  and  mass  transfer  model.  Springer  et  al.  [5] 
predicted  the  net  water-per-proton  flux  ratio  across  the 
membrane.  In  a  later  work.  Springer  et  al.  [6]  fitted  their 
predicted  results  to  the  experimental  data  by  assuming  that 
the  presence  of  liquid  water  only  modifies  the  effective 
porosity  for  oxygen  gas  transport.  Bernardi  and  Verbrugge 
[7,8]  formulated  a  simplified  one-dimensional  model  for 
liquid  water  transport  in  porous  electrodes  assuming  a 
constant  liquid  volume  fraction  and  no  interactions  between 
liquid  and  gas  flows. 

Recently,  Gurau  et  al.  [9]  applied  a  comprehensive  model 
of  fluid  dynamics  and  species  transport  to  a  full  PEM  fuel 
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Nomenclature 

C  species  mass  concentration  (kg  kg-1) 

D  diffusivity  (cm2  s”1) 

F  Faraday  constant  (96,487  C  moP1) 

g  gravitational  acceleration  (cm  s  2) 

hm  mass  transfer  coefficient  between  porous 

cathode  and  gas  channel  (cm  s_1) 

H  total  height  of  cathode  and  channel  (cm) 

Hc  cathode  thickness  (Hc  =  H  —  Hgc)  (cm) 

Hgc  gas  channel  height  (cm) 

I  current  density  (A  cm~2) 

7cr  critical  current  density  (A  cm~2) 

/0  exchange  current  density  (A  cm~2) 

j  mass  flux  in  y-direction  (kg  cm  2  s  ') 

j  mass  flux  vector  (kg  cm-2  s”1) 

krg  relative  permeability  of  gas  phase 

kr\  relative  permeability  of  liquid  phase 

K  permeability  of  porous  cathode  (cm  2) 

L  cathode  length  (cm) 

M  molecular  weight  (kg  moP1) 

p  pressure  (Pa) 

pc  capillary  pressure  (Pa) 

R  gas  constant  (J  moP1  KP1) 

RH  relative  humidity 

s  liquid  water  saturation 

Sh  Sherwood  number  ( hmHgc/D ) 

t  time  (s) 

T  temperature  (°C) 

u  intrinsic  velocity  in  jc-di  recti  on  (cm  s  ') 

u  intrinsic  velocity  vector  (cm  s-1) 

v  intrinsic  velocity  in  y-direction  (cm  s  ') 

Pathode  cathode  potential  versus  SHE  (V) 

Pep  open  circuit  potential  versus  SHE  (V) 

v  coordinate  (cm) 

y  coordinate  (cm) 

Greek  letters 

a  net  water  transport  coefficient  per 

proton 

ac  cathodic  transfer  coefficient 

e  porosity  of  gas  diffusion  cathode 

i]  overpotential  (V) 

p  viscosity  (kg  cnP1  s”1) 

v  kinetic  viscosity  (cm2  s-1) 

9C  contact  angle  (°) 

p  density  (kg  cnP3) 

pk  kinetic  density  (kg  cm-3) 

o  surface  tension  (N  cm-1) 

Superscripts 
H20  water 

02  oxygen 

—  average  value  in  gas  channel 


Subscripts 

flow  flow  channel 

g  gas  phase 

in  inlet 

1  liquid  phase 

ref  reference  value 

sat  saturated  water  vapor 


cell.  Similarly,  Yi  and  Nguyen  [10]  analyzed  two-dimen¬ 
sional  hydrodynamics  and  multicomponent  transport  in  an 
interdigitated  cathode.  It  should  be  noted  that  these  recent 
models  are  valid  only  in  the  absence  of  liquid  water,  and  they 
do  not  account  for  water  condensation  and  evaporation 
phase  change  within  porous  electrodes. 

The  importance  of  two-phase  transport  in  PEM  fuel  cells 
has  been  repeatedly  stressed  in  the  literature  [5,6].  During 
fuel  cell  operation,  especially  at  high  current  densities, 
liquid  water  is  likely  to  appear  in  the  cathode,  resulting 
in  two-phase  transport  phenomena.  The  transport  processes 
then  become  significantly  more  complicated  due  to  the 
coupled  flow  of  liquid  water  and  gaseous  reactants  in  porous 
media  [11,12].  Liquid  water  transport  by  capillary  action, 
dynamic  interaction  between  single-  and  two-phase  zones 
via  evaporation  and  condensation,  and  effects  of  the  phase 
distribution  on  gas  transport  remain  to  be  explored. 

The  present  paper  describes  a  numerical  study  of  gas- 
liquid,  two-phase  flow  and  transport  in  the  air  cathode  of 
PEM  fuel  cells  including  hydrogen  and  direct  methanol  fuel 
cells.  In  the  following  section,  single-  and  two-phase  trans¬ 
port  regimes  are  first  defined  and  discussed  analytically, 
which  is  followed  by  an  accurate  numerical  model  based  on 
the  multiphase  mixture  model  previously  developed  by 
Wang  and  Cheng  [13]  for  two-phase  flow  and  multicompo¬ 
nent  transport  in  porous  media.  Two-dimensional  simula¬ 
tions  are  presented  to  corroborate  the  analytical  estimate  of  a 
threshold  current  density  defining  the  two  transport  regimes 
as  well  as  to  study  the  two-phase  hydrodynamics  and 
transport  of  air  and  liquid  water  and  its  effects  on  the 
performance  of  the  PEM  fuel  cell  cathode. 

2.  Problem  statement 

Consider  a  porous  cathode  adjacent  to  a  flow  channel,  as 
shown  in  Fig.  1.  The  catalyst  layer  is  simplified  as  an 
infinitely  thin  interface  between  the  membrane  and  the 
porous  cathode.  The  oxygen  reduction  reaction  at  the  mem¬ 
brane/cathode  interface  consumes  oxygen  and  generates 
water.  Transport  of  the  reactant  and  product  occurs  in 
the  porous  cathode  and  flow  channel.  As  the  current  density 
(i.e.  the  water  generation  rate)  increases,  water  is  present 
first  as  vapor  only  but  gradually  condenses  to  liquid 
thus  leading  to  the  development  of  two-phase  flow  and 
transport  in  the  air  cathode.  In  the  present  paper,  we  define 
onset  of  the  two-phase  operating  regime  when  liquid 
water  first  appears  in  the  vicinity  of  the  membrane/cathode 


42 


Z.H.  Wang  et  al. /  Journal  of  Power  Sources  94  (2001)  40-50 


interface.  Correspondingly,  a  threshold  current  density  can 
be  identified  for  the  boundary  between  the  single-  and  two- 
phase  transport  regimes  of  the  air  cathode. 

To  obtain  an  analytical  estimate  of  the  threshold  current 
density,  consider  a  one-dimensional  water  transport  process 
along  y-direction  only.  Water  generated  at  the  membrane/ 
cathode  interface  is  transported  through  the  cathode  into  the 
gas  channel  and  then  carried  away  to  the  channel  exit. 
Assuming  that  the  net  water  transport  across  the  membrane 
to  the  cathode  is  quantified  by  a  net  water  transport  coeffi¬ 
cient,  the  rate  of  water  generation  at  the  membrane/cathode 
interface  can  then  be  expressed,  according  to  Faraday’s  law, 
as 


;H2Oi 

J  I  y=H  — 


MH2°(1  +2a)  r 

IF 


(1) 


where  a  is  the  net  water  transport  coefficient  across  the 
membrane,  i.e.  the  ratio  of  the  water  flux  across  the 
membrane  from  the  anode  to  the  proton  flux  transported. 
Physically,  this  water  flux  is  a  combined  result  of  electro- 
osmotic  drag,  fluid  convection,  and  molecular  diffusion. 

In  the  case  of  no  liquid  formation,  gas  diffusion  is 
generally  the  main  mechanism  for  water  transport  through 
the  porous  cathode,  and  gas  convection  is  negligible  due  to 
the  small  permeability.  The  water  vapor  flux  across  the 
cathode  can  thus  be  calculated  by 


H2Oi  _  H2Ol 

;H20 1  _  n  H2Oyg  \y=H  >g  b'=Hgc 

g  I  y=Hgc  -  Ug  S 


H, 


gc 


(2) 


where  p^2° \y=H  represents  the  water  vapor  density  at  the 
membrane/cathode  interface,  and  H  —  Hgc  the  thickness  of 
porous  cathode,  Hc.  The  water  vapor  density  at  the  cathode/ 
channel  interface,  Pg2°  |  v=//  c ,  depends  on  the  mass  exchange 
process  between  the  porous  cathode  and  gas  channel.  Using 
the  average  water  vapor  density  of  the  cross-section  area  in 

the  gas  channel,  p^fl°w,  one  can  express  the  mass  flux  of 
water  vapor  at  the  cathode/channel  interface  as 


■H2Oi 

Jg  \y=Hgc 


hm 


\y=Hgc 


H2°  \ 
Pg,flow  J 


(3) 


where  hm  is  the  mass  transfer  coefficient  between  the  porous 
cathode  and  gas  channel.  Combining  Eqs.  (2)  and  (3)  results 
in 


•H2Oi 

Jg  \y=Hgc 


H20|  _  h2o 

Pg  I  y=H  Pg.flow 


1  /hm+Hc/Df°s 


(4) 


The  two  mass  transfer  resistances  in  the  denominator  of 
Eq.  (4)  are  representative  of  convection  mass  exchange 
between  the  cathode  and  channel  and  diffusion  in  the 
cathode  respectively.  The  water  removal  flux  given  by 
Eq.  (4)  should  equal  the  water  generation  rate  given  by 
Eq.  (1)  under  the  steady-state  condition. 

When  the  water  vapor  density  at  the  membrane/cathode 
interface  reaches  the  saturated  value  corresponding  to  the 
operating  temperature,  liquid  water  begins  to  appear,  mark¬ 
ing  the  end  of  the  single -phase  regime.  Combining  Eqs.  (1) 
and  (4)  thus  yields  the  following  expression  for  the  threshold 
current  density: 


Icr  = 


l/hm  +  tfc/£>g2°£MH2°(l  +  2c 


(5) 


Note  that  the  threshold  current  density  given  by  Eq.  (5)  is 
dependent  on  the  average  water  vapor  density  of  the  cross- 
section  area  in  the  gas  channel.  This  parameter  increases 
down  the  channel  due  to  the  continual  addition  of  water 
vapor  from  the  cathode;  namely, 


fin  //gc 


Pg.flow 


MH2°(  1  +2a) 


(6) 


With  Eq.  (6),  Eq.  (5)  can  now  be  rewritten  as 

2 F  Pg,2s°(l  -  RHin)  (  L  11  Hc\~l 
Cl  MTO  1  +  2a  \Hgcuin  hm  Df°e) 

(V) 


where 


RHin  = 


pH2° 

/_gan_ 

H2° 

r  g,sat 


(8) 


is  the  relative  humidity  of  gas  at  the  inlet.  Note  that  the 
threshold  current  density  is  dependent  on  the  water  vapor 
diffusion  resistance  in  the  porous  cathode,  the  mass  transfer 
resistance  at  the  porous  cathode/gas  channel  interface,  and 
the  ability  of  air  in  the  channel  to  carry  away  water  vapor.  It 
reduces  with  the  increases  of  inlet  relative  humidity,  channel 
length,  and  porous  cathode  thickness  and  the  decreases  of 
inlet  velocity,  operating  temperature,  and  channel  height. 

All  parameters  required  to  estimate  the  threshold  current 
density  by  Eq.  (7)  are  readily  available  from  the  operating 
conditions  except  for  the  net  water  transport  coefficient  and 
the  mass  transfer  coefficient  at  the  interface  between  the 
porous  cathode  and  flow  channel.  The  net  water  transport 
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Table  1 


Property  parameters  used  in  numerical  simulations 


Parameter 

Symbol 

Value 

Reference 

Water  vapor  diffusivity 

DWO 

0.355  cm2  s-1 

Cussler  [15] 

Oxygen  diffusivity 

Dk 

0.284  cm2  s~ 1 

Cussler  [15] 

Porosity  of  air  cathode 

£ 

0.3 

Yi  and  Nguyen  [10] 

Net  water  transport  coefficient 

a 

0.3 

Assumed 

Permeability  of  porous  cathode 

K 

1.0  x  10~7  cm2 

Assumed 

Surface  tension 

G 

6.25  N  cm-1 

Incropera  and  DeWitt  [16] 

Constant  in  Tafel  equation 

Vcg°jef 

0.3125  A  cm^2 

Yi  and  Nguyen  [10] 

Cathodic  transfer  coefficient 

ac 

2 

- 

Open  circuit  potential 

^ocp 

1.1  V 

Yi  and  Nguyen  [10] 

coefficient  is  a  combined  result  of  electro-osmotic  drag  by 
electric  field,  fluid  convection  by  pressure  difference,  and 
molecular  diffusion  by  species  concentration  difference  in 
the  membrane  region,  thus  its  rigorous  determination 
requires  a  full  cell  analysis.  Some  estimate  of  this  parameter 
has  been  made  by  Yi  and  Nguyen  [10]  for  hydrogen  fuel 
cells. 

To  estimate  the  mass  transfer  coefficient  at  the  cathode/ 
channel  interface,  consider  mass  transfer  in  a  fully-devel¬ 
oped  laminar  flow  through  a  parallel-plate  channel  with  a 
constant  mass  flux  applied  at  one  surface  and  no-flux  at  the 
other.  In  this  situation,  the  dimensionless  mass  transfer 
coefficient  defined  by  Sherwood  number  is  given  by  [14] 

Sh  =  =  2.693  (9) 

Dg 

Although,  the  fuel-cell  flow  channel  differs  slightly  from  the 
geometry  described  above  because  the  cathode/channel 
interface  is  permeable,  the  calculation  based  on  the  detailed 
numerical  results  to  be  reported  in  the  following  yields  a 
Sherwood  number  of  2.637  for  the  present  problem.  This 
good  agreement  with  Eq.  (9)  not  only  indicates  that  there  is 
negligibly  small  flow  through  the  cathode  of  low  perme¬ 
ability,  but  also  confirms  that  the  interfacial  mass  transfer 
coefficient  can  be  safely  represented  by  Eq.  (9). 

Using  the  properties  listed  in  Table  1  and  the  operating 
parameters  summarized  in  Table  2,  Fig.  2  illustrates  the 
effect  of  the  air  inlet  velocity  on  the  threshold  current 
density.  It  can  be  seen  that  the  threshold  current  density 

Table  2 


Base  case  and  its  operating  conditions 


Parameter 

Symbol 

Value 

Total  height  of  cathode  and  channel 

H 

0.12  cm 

Channel  height 

tfgc 

0.07  cm 

Cathode  length 

L 

2cm 

Relative  humidity  at  the  inlet 

RHln 

3.43% 

Oxygen  concentration  at  the  inlet 

r^Oi 

g,in 

0.21 

Air  velocity  at  the  inlet 

Um 

30cm  s-1 

Operating  pressure 

P 

latm 

Operating  temperature 

T 

80°C 

Cathode  overpotential 

n 

-0.3473  V 

approaches  an  asymptotic  value  at  um  of  about  2ms'  , 
implying  that  the  water  vapor  removal  is  limited  by  diffusion 
through  the  porous  cathode  under  this  condition. 

While  the  above  approximate  analysis  and  classification 
yield  simple  quantification  of  the  two  regimes  of  water 
transport  and  distribution  in  the  air  cathode,  accurate  mod¬ 
eling  of  the  two-phase  flow  and  transport  must  be  performed 
numerically  as  is  shown  in  Section  3. 

3.  Numerical  modeling 

Several  prior  models  have  treated  the  single-phase  opera¬ 
tion  of  the  air  cathode  [9,10,17].  The  two-phase  regime 
involving  two-phase  flow  and  transport  in  the  air  cathode  has 
yet  to  be  modeled.  Moreover,  no  model  has  been  able  to 
include  both  the  single-  and  two-phase  regimes  as  well  as 
their  transition.  Therefore,  the  main  objective  of  this  section 
is  to  present  a  unified  model  that  encompasses  both 
the  single-  and  two-phase  regimes,  and  ensures  a  smooth 


Fig.  2.  Effect  of  the  air  inlet  velocity  on  the  threshold  current  density 
defining  the  boundary  between  the  single-  and  two-phase  regimes. 
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transition  between  the  two.  To  account  for  liquid  phase 
transport  by  capillary  action  in  the  two-phase  zone,  evapora¬ 
tion-condensation  phase  change  at  the  evaporation  front, 
and  multicomponent  gas  transport  due  to  both  diffusion  and 
convection  within  the  single -phase  region,  the  multiphase 
mixture  model  previously  developed  by  Wang  and  co-work- 
ers  [13,18-22]  is  conveniently  applied.  The  reader  is 
referred  to  Wang  and  Cheng  [13]  for  details  of  the  multi¬ 
phase  mixture  model  and  its  applications  to  a  number  of 
multiphase  transport  problems  in  porous  media.  In  the 
following,  only  the  model  assumptions  and  governing 
equations  for  the  present  problem  are  summarized. 


Here,  the  momentum  balance  equation  has  been  modified  to 
be  valid  in  both  the  porous  cathode  and  the  open  flow 
channel,  reducing  to  the  extended  Darcy’s  law  for  two-phase 
flow  in  porous  media  with  a  small  permeability  and  the 
Navier-Stokes  equation  inside  the  flow  channel  with  the 
porosity  being  unity  and  the  permeability  being  infiniteness. 
Note  that  u  is  the  intrinsic  velocity  vector  based  on  the  open 
pore  area  only.  The  species  transport  equation  is  applicable 
to  oxygen  (02)  as  well  as  water  (H20).  Also,  the  species 
concentration  on  the  mass  basis  has  been  used  instead  of  the 
partial  density  employed  in  the  foregoing  order-of-magni- 
tude  analysis.  They  are  related  by 


3.1.  Model  assumptions 

A  real  gas-diffusion  electrode  has  an  extremely  complex 
geometric  structure,  featuring  both  micro-  and  macro-pores 
as  well  as  a  distribution  of  pore  sizes.  In  addition,  the  liquid 
water  distribution  and  the  effective  surface  tension  of  the 
gas-liquid  interface  within  a  solid  matrix  depend  on  hydro- 
phobic  or  hydrophilic  nature  of  the  electrode.  As  a  first  step 
towards  modeling  the  two-phase  distribution  and  transport 
in  air  cathodes,  it  is  necessary  to  invoke  simplifying  assump¬ 
tions  that  are  commonly  made  in  the  current  literature.  These 
are  as  follows. 

1.  The  gas  diffusion  electrode  is  isotropic  and  homoge¬ 
neous,  and  characterized  by  an  effective  porosity  and 
permeability. 

2.  The  catalyst  layer  is  treated  as  an  infinitely  thin  surface 
of  reaction. 

3.  The  cell  temperature  remains  constant. 

4.  The  gas  phase  is  an  ideal  mixture. 

In  future  work,  efforts  will  be  made  to  incorporate  more 
realistic  characteristics  of  the  air  cathode  structure  into  the 
present  model  framework. 

3.2.  Governing  equations 
Continuity 

^l+\7  .(epu)  =  0  (10) 

Momentum  conservation 

d(epu)  7  p 

— - - h  V  •  ( spuu )  =  V  •  (epvw)  —  evp  +  spkg  —  a  —u 

at  K 

(11) 


Species  conservation 
d 

—  (epC)  +  V  •  ( eycpuC ) 

=  V  •  (epDVC)  +  V  •  (efpj&DiVCi 
+  pg(  1  -  s)DgVCg  -  pDVC]}  -  V  •  [(Q  -  Cg»] 

(12) 


P'i 

C\  =  -1  (13) 

7  Pj 

P1 

C'=-  (14) 

P 

where  superscript  i  represents  a  species  of  02  or  H20,  and 
subscript  j  a  gas  (g)  or  liquid  (1)  phase. 

A  unique  feature  of  the  multiphase  mixture  model  is  that  it 
does  not  need  to  track  phase  interfaces  separating  single  - 
from  two-phase  regions,  and  hence  greatly  simplifies  nume¬ 
rical  simulation  of  the  present  problem.  In  fact,  all  the  above 
governing  equations  identically  reduce  to  their  single-phase 
counterparts  in  the  limits  of  the  liquid  saturation,  s,  equal  to 
zero  and  unity,  respectively. 


3.3.  Mixture  parameters 


In  the  governing  Eqs.  ( 10)— (12),  the  mixture  variables  and 
properties  are  defined  as 

Density:  p  =  pxs  +  pg(l  —  5)  (15) 

Concentration  :  pC  =  p\C\S  +  pgCg(l  —  s)  (16) 

Velocity  :  pu  =  pxU\  +  p„ug  (17) 

Kinetic  density  :  Pk  =  PiM(s)  +  pg2g(y)  (18) 

PiS  +  pJl-s) 


Viscosity  : 


P  = 


(19) 


(kvl/vi)  +  (£rg/vg) 

Diffusion  coefficient  :  pD  =  pxsD\  +  p  (1  —  s)D„ 

(20) 

p{hC\  +  AgCg) 


Advection  correction  factor  : 


7c  = 


Relative  mobilities  : 


2i(5)  = 


P\sC\  +  pg(l  -  s)Cg 

(21) 

Wvi 


WV1  +  £rg/vg  ’ 

Ag(s)  =  1 -^(j)  (22) 

In  addition  to  the  above  conservation  equations,  the  multi¬ 
phase  mixture  formulation  embodies  the  following  impor¬ 
tant  relation  for  calculating  the  individual  phase  velocities 
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from  the  mixture  flow  field: 

j\  =  [V/?c  +  (Pi  -  Pg)g]  (23) 

Physically,  the  first  term  on  RHS  represents  the  capillary 
flow  due  to  the  gradient  in  meniscus  curvature,  whereas  the 
second  term  accounts  for  gravity-induced  phase  migration. 
The  individual  phase  velocities  (i.e.  intrinsic  velocities 
based  on  the  open  pore  area  rather  than  phase-occupied 
areas)  are  then  calculated  from 

epjMi  =  j\  +  2.\epu  (24) 

epgMg  =  — y'|  +  Agspu  (25) 


be  obtained  from  the  steam  table.  Thus,  the  liquid  saturation, 
s,  can  be  determined  from  the  mixture  concentration  of 
water,  C,  via  the  following  relation: 


s 


p,(c ; 


h2o 


p8(c*°  -  egg) 
-  cH=°)  +  pg(CH=° 


rH2ON 

'-'g,sat/ 


(29) 


Similarly,  for  oxygen 


ro2  _ 

l-'l,sat 


0  Cg,t 


Pig 

Pg(l  —  S) 


+  1 


C°2 


(30) 


within  the  two-phase  zone,  assuming  that  oxygen  is  inso¬ 
luble  in  the  liquid  phase. 


3.4.  Constitutive  relations 

As  a  first  approximation,  the  gas-diffusion  electrode  is 
treated  as  a  homogeneous  porous  medium  with  the  relative 
permeabilities  for  liquid  and  gas  phases  represented  by 
[18] 

k{\  =  s 3  and  kIg  =  (1  —  5) 3  (26) 

Similarly,  the  capillary  pressure  between  the  two  phases  can 
be  expressed  as  [13] 

\  1/2  , 

J  [1.417(1  —  s)  —  2.120(1  —  s)~ 

+  1.263(1  —  v)3]  (27) 

where  er  is  the  surface  tension  at  the  gas-liquid  interface,  and 
9C  the  contact  angle.  The  surface  tension  is  taken  to  be 
0.0625  Nm~’  for  the  liquid  water-air  system  at  80°C, 
whereas  9C  is  dependent  upon  the  hydrophilic  (i.e. 
0°  <  9C  <  90°)  or  hydrophobic  (i.e.  90°  <  9C  <  180°)  nat¬ 
ure  of  the  gas-diffusion  layer,  and  thus  varies  with  the  Teflon 
content.  In  the  present  work,  9C  is  chosen  to  be  zero  to 
simulate  a  hydrophilic  electrode.  A  parametric  study  by 
varying  9C  will  be  conducted  in  a  separate  publication  to 
access  important  effects  of  the  wetting  property  of  air 
cathodes  on  two-phase  transport  and  distribution. 


Pc  =  0  cos  9C  y- 


3.6.  Boundary  conditions 


At  the  inlet  of  the  gas  channel,  constant  flow  rate  and  mass 
concentration  of  each  species  are  specified.  At  the  outlet, 
both  velocity  and  concentration  fields  are  assumed  to  be 
fully  developed.  Hence,  at  0  <  y  <  H„c 


U\x=0  ~  VL=0  —  O’ 

/^02 1  _  W32 

C'  lx=0  —  Hn 

=  V\x=L  = 

=  0 


du 

dx\ x=L 

dC° 


dx  \x=L 

at  H„c  <  y  <  H 

u\x=0  =  VL=0  =  °» 

dC° 2 


dx 


=  0 


x=0 


U\ x=L  =  V\x=L  = 

dC(h 


dx 


=  0 


C 


H2Op  TT 

.H2Oi  __  /'g.sat^nir 


\x=0 


dc^-o 


dx 


=  0, 


x=L 


dC 


H,0 


dx 


dCH2 


=  0, 


x=0 


dx 


=  0, 


x=L 


(31) 


(32) 


(33) 


(34) 


x =L 


3.5.  Equilibrium  conditions 


At  the  membrane/cathode  interface  where  the  oxygen  reduc¬ 
tion  reaction  occurs 


In  the  present  study,  liquid  water  is  considered  to  appear 
in  the  porous  cathode  only  when  the  water  vapor  pressure 
reaches  its  saturated  value  at  the  operating  temperature. 
Inside  the  two-phase  zone,  thermodynamic  equilibrium  is 
assumed  to  hold  true,  and  thus  the  concentrations  of  water  in 
gas  and  liquid  phases  given  by  their  equilibrium  values, 
respectively,  are 


rH2o 

^g,sat 


MH2V(r) 

RpaT 


Cg?  =  l 


(28) 


where  pv(T)  is  the  water  vapor  saturation  pressure,  and  can 


U\ y=H  ~  0; 

and 


y=H  = 


(J02  +/hO)| 


y=H 


P\y=H 


(35) 


ycpsvC-r  y  +  (^  -  Cg.sat)  (p\-pg)g 

=  j\y=H  (36) 

for  both  species  O2  and  H20.  Eq.  (36)  was  derived  by  Cheng 
and  Wang  [22]  to  include  three  diffusion  mechanisms  in  a 
compact  form:  mass  diffusion  through  liquid  and  gas  phases, 
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respectively,  and  capillary  diffusion  within  the  two-phase 
zone.  The  effective  diffusivity,  T,  in  Eq.  (36)  is  thus  defined 
as 


r  =  aDg  in  all  gas  regions, 

r  =  (  L^t - SMt)Pi — y  ;n  the  two-phase  region  (37) 

Pi  6-  l,sat  Pg^g,sat 


where  the  capillai'y  diffusion  coefficient,  Dc,  is  given  by 


(  _  jdjgg  (  dpc\ 

v  \  ds  J 


(38) 


The  02  mass  flux  at  this  boundai'y  can  be  expressed  as 
M°2 


i° 2  I 

J  I  y=H 


4F 


(39) 


whereas  the  EEO  mass  flux  has  been  given  by  Eq.  (1).  Tafel 
equation  is  employed  to  describe  the  local  current  density 
along  the  membrane/cathode  interface;  namely, 


/  =  /„ -  exp 

Cg,ref 


(40) 


where  the  term  (1  —  s)  is  used  to  account  for  the  fraction  of 
surface  rendered  inactive  by  the  presence  of  liquid  water. 
At  the  bottom  of  the  gas  channel, 


U\y=0  =  0. 

dC° 2 


b=o 


=  0, 


<9CH2° 


dy 


=  0, 


>’=0 


dy 


=  0 


(41) 


y=0 


3.7.  Cathode  potential 

Based  on  the  reference  that  the  membrane  phase  potential 
is  zero  at  the  membrane/cathode  interface,  the  cathode 
potential  can  then  be  related  to  the  cathode  polarization 
oveipotential  by 

^cathode  —  C^ocp  4”  t)  (42) 

where  V(,cp  is  the  open  circuit  potential  of  the  air  cathode  in 
reference  to  the  standard  hydrogen  electrode  (SHE).  Note 
that  the  cell  voltage  is  equal  to  the  cathode  potential  minus 
the  ohmic  potential  drop  across  the  membrane  and  the  anode 
oveipotential. 


4.  Results  and  discussion 

The  coupled  governing  Eqs.  (10)— (12)  with  the  boundary 
Eqs.  (3 1)— (36)  and  (41)  are  solved  numerically  by  a  finite- 
volume-based  finite  difference  method.  Numerical  details 
were  previously  given  by  Wang  [18],  and  thus  are  not 
repeated  here.  The  relevant  parameters  used  in  all  simula¬ 
tions  reported  in  the  following  are  listed  in  Table  1. 


Fig.  3.  Polarization  curve  encompassing  both  single-  and  two-phase 
regimes  of  the  porous  cathode. 


Fig.  3  shows  a  simulated  cathode  polarization  curve  under 
the  operating  conditions  described  in  Table  2.  Note  that  a  dry 
inlet  is  chosen  as  a  case  of  study  here  because  of  recent 
interest  in  PEM  fuel  cells  without  external  humidification 
[23]  or  low-humidity,  self-sustaining  fuel  cells  for  portable 
electronics.  In  addition,  the  dry  inlet  of  the  air  cathode  is  of 
potential  relevance  to  direct  methanol  fuel  cells  [24]. 

The  cathode  polarization  curve  is  obtained  by  varying  the 
surface  oveipotential  with  a  constant  air  flow  rate  at  the  inlet. 
This  curve  only  considers  the  voltage  loss  at  the  cathode,  and 
does  not  include  the  anode  kinetic  polarization  and  mem¬ 
brane  ohmic  loss.  Thus,  the  polarization  curve  is  higher  than 
experimental  results  as  expected,  but  its  general  trend 
compares  well  to  typical  experiments  [6].  At  low  current 
densities,  the  cathode  activation  oveipotential  is  solely 
responsible  for  the  potential  losses  of  the  cathode.  At  higher 
current  densities,  more  oxygen  is  consumed  and  more  water 
is  generated  at  the  membrane/cathode  interface  due  to  both 
the  electrochemical  reaction  and  the  water  transport  across 
the  membrane  from  the  anode.  The  cathode  polarization 
then  begins  to  be  limited  by  the  transport  of  02.  As  the 
current  density  exceeds  the  threshold  value,  i.e.  1.35  A  cm~2 
in  the  present  case  (in  comparison,  its  analytical  estimate  is 
1.26  A  cm~2  as  can  be  obtained  from  Fig.  2),  more  water  is 
produced  than  that  possibly  carried  away  by  air  from  the 
electrode,  and  liquid  water  starts  to  form  in  the  vicinity  of  the 
membrane/cathode  interface.  The  predicted  onset  of  the 
two-phase  regime  at  1.35  A  cm  2  may  appear  a  little  high, 
but  this  value  corresponds  to  an  extremely  small  cell  (i.e. 
2  cm  long)  and  the  dry  inlet  air.  Both  Eq.  (7)  and  numerical 
calculations  indicate  that  the  onset  current  density  for  the 
two-phase  regime  drops  to  0.64  A  cm  2  for  5  cm  long  cells 
with  dry  air  at  the  inlet. 

The  presence  of  liquid  water  reduces  the  pore  spaces  for 
oxygen  transport  to  the  reaction  surface  as  well  as  renders 
part  of  the  surface  electrochemically  inactive.  These  effects 
result  in  a  slightly  steeper  slope  of  the  cathode  potential 
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curve  upon  the  threshold  current  density,  as  can  be  seen  from 
Fig.  3.  Clearly,  a  two-phase  model  is  warranted  for  fuel  cells 
operated  at  such  high  current  densities,  an  operating  range  of 
special  interest  for  transportation  applications.  An  advan¬ 
tage  of  the  present  two-phase  model  is  that  it  not  only 
provides  a  unified  formulation  for  both  single-phase  analysis 
at  low  current  densities  and  two-phase  calculation  at  higher 
values,  but  also  predicts  and  ensures  the  smooth  transition 
from  the  single-  to  two-phase  region  automatically. 

With  further  increase  in  the  cell  current  density,  the  two- 
phase  zone  expands  and  its  evaporation  front  propagates 
towards  the  cathode/channel  interface.  The  numerical  simu¬ 
lation  shows  that  the  evaporation  front  reaches  the  cathode/ 
channel  interface  at  an  average  current  density  of 
1 .47  A  cm  2,  beyond  which  point  two-phase  flow  appears 
in  the  cathode  channel.  This  is  beyond  the  scope  of  the 
present  work. 

Fig.  4  shows  the  velocity  field  of  the  two-phase  mixture  in 
the  cathode  and  gas  channel  for  the  base  case  (i.e.  at  the 
average  cell  current  density  of  1.4  A  cm  2).  As  expected, 
there  is  a  large  difference  in  the  velocity  scale  between  the 
porous  region  and  the  open  channel.  The  mixture  velocity  in 
the  porous  cathode  is  at  least  two  orders  of  magnitude 
smaller  than  that  in  the  open  channel,  indicating  that  gas 
diffusion  is  the  dominant  transport  mechanism  in  the  single¬ 
phase  region  of  the  porous  cathode.  The  flow  field  in  the 
open  channel  is  fully  developed  in  view  of  the  large  aspect 
ratio  of  the  channel  length  to  width  (equal  to  29  in  the 
present  case),  as  can  be  seen  from  Fig.  4  where  the  channel 
length  is,  however,  not  drawn  to  scale  for  better  view.  The 
velocity  distribution  in  the  conventional  flow  field  such  as 
the  one  studied  in  the  present  work  is  relatively  simple  and 
could  be  analytically  determined;  however,  a  numerical 
analysis  based  on  the  comprehensive  model  as  presented 
here  becomes  necessary  for  complex  designs  such  as  the 
interdigitated  flow  field  [10].  Magnified  plots  of  individual 


1.5m/s 


Fig.  4.  Two-phase  mixture  velocity  field  in  the  porous  cathode  and  flow 
channel  at  the  current  density  of  1.4  A  cm~2. 


Fig.  5.  Water  vapor  concentration  contours  in  the  porous  cathode  and  flow 
channel  at  the  current  density  of  1.4  A  cm~2. 


phase  velocity  fields  within  the  porous  cathode  region  are 
presented  in  Fig.  7(a)  and  (b)  and  discussed  in  detail  therein. 

Fig.  5  displays  the  water  vapor  concentration  distribution 
in  the  porous  cathode  and  flow  channel.  As  the  air  flows 
down  the  channel,  water  vapor  is  continually  added  from  the 
cathode,  resulting  in  an  increased  water  vapor  concentration 
along  the  channel.  This  leads  to  a  decreased  concentration 
gradient  in  water  vapor,  and  hence  a  lower  water  vapor 
diffusion  flux  from  the  membrane/cathode  surface  to  the  gas 
channel.  As  a  result,  liquid  water  may  first  appear  in  the 
vicinity  of  the  membrane/cathode  interface  near  the  channel 
outlet.  A  two-phase  zone  at  this  location  is  indeed  predicted 
in  the  present  simulation  shown  in  Fig.  5,  where  the  water 
vapor  mass  concentration  is  seen  to  be  kept  at  the  saturated 
value,  i.e.  0.2907,  corresponding  to  the  operating  tempera¬ 
ture  of  80°C.  It  is  also  seen  that  the  two-phase  zone  is 
confined  within  the  cathode  at  the  average  current  density  of 
1.4  A  cm~2. 

Fig.  6  shows  the  liquid  water  saturation  distribution  in  the 
same  case.  In  accordance  with  Fig.  5,  liquid  water  is  seen  in 
the  upper-right  comer  to  coexist  with  the  saturated  water 
vapor.  The  largest  liquid  amount  predicted  is  around  6.3%  at 
the  current  density  of  1.4  A  cm  2.  This  value  is  lower  than 
that  reported  experimentally,  likely  because  the  present 
study  simulates  a  single  channel  whereas  experimental  fuel 
cell  fixtures  usually  use  series-connected  channels.  With  the 
effective  channel  length  extended,  a  higher  liquid  saturation 
may  result.  Also,  a  dry  inlet  is  used  in  the  present  simula¬ 
tions,  whereas  many  experimental  systems  humidify  the  air 
stream.  The  numerical  results  indicate  that  the  mass  flux  of 
liquid  water  due  to  gravity  is  less  than  0. 1  %  of  that  caused  by 
capillary  action  in  the  two-phase  zone.  This  is  because  the 
Bond  number  [16],  defined  as  (px  —  p„)gH^ /a,  is  about 
0.04,  implying  the  negligible  gravity  effect  as  compared 
with  the  surface  tension  effect  in  the  two-phase  zone.  Thus, 
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Fig.  6.  Liquid  water  saturation  contours  in  the  porous  cathode  and  flow 
channel  at  the  current  density  of  1.4Acm~2.  The  evaporation  front 
separating  the  two-phase  zone  from  the  single-phase  region  is  approxi¬ 
mately  represented  by  the  contour  of  s  =  0.01. 

the  liquid  water  is  transported  mainly  by  capillary  action  to 
the  evaporation  front.  The  front  can  be  approximately 
depicted  by  the  contour  of  s  =  0.01  shown  in  Fig.  6. 


To  display  individual  velocity  field  of  each  phase  in  the 
porous  cathode,  the  gas  and  liquid  velocities  are  calculated 
from  the  mixture  velocity  according  to  Eqs.  (24)  and  (25) 
and  shown  in  Fig.  7(a)  and  (b),  respectively.  Within  the  two- 
phase  zone,  liquid  water  flows  towards  the  evaporation  front 
under  capillary  action.  The  maximum  liquid  velocity  is 
approximately  2.3  x  10“5  m  s_1.  The  gas  velocity  field 
shown  in  Fig.  7(b)  differs  from  the  liquid  velocity  by  three 
orders  of  magnitude  and  displays  opposite  patterns  in  the 
single-  and  two-phase  regions.  In  the  single -phase  region, 
the  gas  velocity  is  primarily  normal  towards  the  cathode/ 
channel  interface  due  to  the  non-zero  gas  velocity  created  by 
the  electrochemical  reaction  of  H+/02.  It  is  noted  from 
Eq.  (35)  along  with  Eqs.  (1)  and  (39)  that  the  normal  velocity 
at  the  reaction  surface  is  negative  (pointing  away  from  it),  in 
accordance  with  the  physical  process  that  hydrogen  ions  are 
recombined  with  02  to  produce  water.  Across  the  evapora¬ 
tion  front,  however,  the  gas  flow  reverses  its  direction, 
leaving  away  from  the  evaporation  front  again  due  to 
capillary  action  in  the  two-phase  zone.  The  maximum  gas 
velocity  is  approximately  0.005  m  s~ 1  in  the  single-phase 
region,  while  0.02  ms-1  in  the  two-phase  zone.  For  a 
cathode  thickness  of  0.0005  m  and  the  gas  diffusivity  of 
around  3  x  10“5  m2  s  1  for  both  water  vapor  and  oxygen, 
the  Peclet  number  for  gas  transport  is  thus  0.08  and  0.32  for 
single-  and  two-phase  zones,  respectively.  This  implies  a 
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Fig.  7.  Liquid  and  gas  velocity  fields  at  the  average  current  density  of  1.4  A  cm  2:  (a)  liquid;  (b)  gas. 
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Fig.  8.  Oxygen  concentration  contours  in  the  porous  cathode  and  flow 
channel  at  the  current  density  of  1.4  A  cm~2. 

definitely  negligible  contribution  of  gas  convection  in  the 
single -phase  region  but  some  effect  of  gas  convection  in  the 
two-phase  zone. 

Fig.  8  shows  the  oxygen  concentration  distribution  in  both 
the  cathode  and  open  channel  at  the  same  current  density  of 
1.4Acm~2.  It  can  be  seen  that  about  40%  of  the  inlet 
oxygen  has  been  consumed  at  the  half-length  of  the  channel. 
At  higher  current  densities,  this  fraction  becomes  larger, 
leaving  the  channel  outlet  section  nearly  depleted  of  oxygen. 
A  direct  consequence  of  the  non-uniform  oxygen  concen¬ 
tration  profile  along  the  reaction  surface  is  the  non-uniform 
current  density  distribution,  because  the  oxygen  reduction 
reaction  is  first-order  as  indicated  by  Eq.  (40).  This  is  shown 
in  Fig.  9,  where  the  local  current  density  is  equal  to 
2.22  A  cnrf2  at  the  inlet,  but  only  0.86  A  cnrf2  near  the 


outlet.  Thus,  the  water  generated  at  the  inlet  is  2.5  times  of 
that  at  the  outlet.  A  slightly  abrupt  change  in  the  slope  is  also 
noted  from  Fig.  9  at  x  =  0.8  cm  for  the  local  current  density 
profile.  This  location  corresponds  exactly  to  the  beginning  of 
the  two-phase  zone  along  the  reaction  surface.  The  local 
current  density  is  reduced  as  liquid  water  forms,  and  sub¬ 
sequently  impedes  the  oxygen  supply  to  the  electrochemical 
reaction  surface. 


5.  Conclusions 

Single-  and  two-phase  regimes  of  water  transport  and 
distribution  in  PEMFC  air  cathodes  have  been  defined  by  a 
threshold  current  density  corresponding  to  first  appearance 
of  liquid  water  at  the  membrane/cathode  interface.  An 
analytical  estimate  of  the  threshold  current  density  was 
obtained.  A  two-phase  flow  and  transport  model  was  devel¬ 
oped  to  predict  liquid  water  formation  and  its  effects  on  the 
electrochemical  kinetics  at  the  membrane/cathode  interface 
and  reactant  and  product  transport  in  the  air  cathode.  The 
present  model  is  novel  in  that  it  encompasses  both  single  - 
and  two-phase  regimes  corresponding  to  low  and  high 
current  densities  and  therefore  is  capable  of  predicting 
the  transition  between  the  two  regimes.  The  water  liquid 
and  vapor  transport  is  controlled  by  capillary  action  and 
molecular  diffusion,  respectively,  due  to  negligible  small  air 
velocity  within  the  porous  electrode. 

The  present  two-phase  model  is  directly  applicable  to  the 
air  cathode  of  both  hydrogen  and  direct-methanol  fuel  cells. 
In  the  latter  case  (DMFC),  the  two-phase  region  might  be 
more  extensive  due  to  large  electro-osmotic  and  diffusion 
fluxes  of  liquid  water  from  the  DMFC  anode.  Similarly,  the 
DMFC  anode,  where  the  liquid  phase  is  a  methanol  aqueous 
mixture  and  the  gas  phase  is  carbon  dioxide,  can  be  treated 
by  the  present  model. 

Future  work  will  include  the  following  activities. 

1.  Integrate  the  two-phase  cathode  model  into  a  multi¬ 
dimensional  full-cell  model  (Um  et  al.  [17])  to  examine 
overall  cell  performance  at  high  current  densities.  This 
integration  will  permit  addressing  the  coupling  between 
the  water  content  and  water  activity  determined  by  the 
state  of  water  at  the  membrane/cathode  interface.  In 
other  words,  the  liquid  water  formation  within  the 
cathode  may  in  turn  impact  the  water  transport  within 
the  membrane  and  hence  its  ionic  conductivity  and 
ohmic  potential  drop. 

2.  Use  the  two-phase  model  to  properly  assess  the 
beneficial  effects  of  an  interdigitated  flow  field  design 
on  cell  performance  in  the  high  current  density  regime. 

3.  Perform  quantitative  comparisons  between  the  numer¬ 
ical  predictions  and  localized  experimental  measure¬ 
ments. 

4.  Extend  the  two-phase  model  to  more  accurately 
incorporate  wetting  characteristics  of  the  gas-diffusion 
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electrode  and  hence  assess  their  effects  on  cell 
performance  at  high  current  densities. 
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